arXiv: 1503.06569v2 [math.NA] 21 May 2015 


NUMERICAL EVALUATION OF TWO AND THREE PARAMETER 
MITTAG LEFFLER FUNCTIONS* 

ROBERTO GARRAPPAt 

Abstract. The Mittag-Lefiler (ML) function plays a fundamental role in fractional calculus but 
very few methods are available for its numerical evaluation. In this work we present a method for the 
efficient computation of the ML function based on the numerical inversion of its Laplace transform 
(LT): an optimal parabolic contour is selected on the basis of the distance and the strength of 
the singularities of the LT, with the aim of minimizing the computational effort and reduce the 
propagation of errors. Numerical experiments are presented to show accuracy and efficiency of the 
proposed approach. The application to the three parameter ML (also known as Prabhakar) function 
is also presented. 
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1. Introduction. The Mittag-Leffler (ML) function was introduced, at the be¬ 
ginning of the twentieth century, by the Swedish mathematician Magnus Gustaf 
Mittag-Leffler [24, 25] while studying summation of divergent series. Extensions 
to two [41] and three [31] parameters of the original single parameter function were 
successively considered; all these functions can be regarded as special instances of the 
generalized hypergeometric function investigated by Fox [9] and Wright [43]. 

Until the 1960s, few authors (e.g., [19]) recognized the importance of the ML 
function in fractional calculus and, in particular, in describing anomalous processes 
with hereditary effects [1, 4, 5, 8]. For an historical outline and a review of the main 
properties of the ML function we refer to [17, 23] and to the recent monograph [15]. 

For any argument z G C, the ML function with two parameters a, /3 G C, with 
5ft(a) > 0, is defined by means of the series expansion 


( 1 . 1 ) 




where r(z) = is the Euler’s gamma function; since the integral repre¬ 

sentation of r(z) holds only for 3ft(z) > 0, the extension to the half-plane 3ft(z) < 0, 
with z ^ {0, —1, —2 ,...}, is accomplished by means of the relationship r(z + n) = 
z{z+!)■■■ {z + n- l)r(z), n e N, [20, 22]. 

As a special case, the ML function with one parameter is obtained for /3 = 1, i.e. 
Eol{z) = Ea,i{z), whilst the generalization to a third parameter 7 


( 1 . 2 ) 


KA^) 


1 ^ r(7 + fc)z*-- 

^A:!r(afc-b/3) 


is recently receiving an increasing attention due to the applications in modeling po¬ 
larization processes in anomalous or inhomogeneous materials [3]. 

In this work we restrict our attention to real parameters a, P and 7, with a > 0 
and 7 > 0, since they are of more practical interest. 
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With the exception of a few special cases in which the ML function can be repre¬ 
sented in terms of other elementary and special functions, as for instance Ei^i[z) = e^, 
E 2 ,i{z^) = cosh( 2 :), i? 2 ,i(—2^) = cos( 2 :) and ^(± 2 :^/^) = e^erfc(=F-2^^^), most of the 
programming languages do not provide built-in functions for the ML function. 

Although it is theoretically possible to evaluate Ea,p{z) and p{z) by truncating 
the series in (1.1) and (1.2), in the majority of the applications this is not advisable: 
with arguments having moderate or large modulus \z\ the convergence of the series 
is very slow and involves an exceedingly large amount of computation; moreover, a 
large number of terms in the series can significantly grow before decreasing, thus 
generating overflow or numerical cancellation unless variable precision arithmetic is 
used. In finite precision arithmetic (which is the natural environment for scientific 
computing) the use of (1.1) and (1.2) is, therefore, confined just to small arguments. 

Very few methods have been presented so far in the literature. The sophisticated 
algorithm described in [16] uses different techniques to evaluate the ML function and 
its derivative in different parts of the complex plane. Other approaches based on 
mixed techniques (Taylor series, asymptotic series, and integral representations) were 
discussed in [18, 32]. The only existing Matlab code [30] (which implements some of 
the ideas introduced in [16]) shows a great variability in the amount of computation 
required to achieve a prescribed accuracy and in some regions of the complex plane 
turns out to be poorly accurate. 

The recent introduction [11, 12, 26] of new methods for fractional differential 
equations involving a large number of evaluations of the ML function, also with matrix 
arguments [27], motivates the investigation of different techniques to perform the 
computation, over the whole complex plane, in an accurate and fast way. 

In this paper we consider an approach based on the inversion of the Laplace 
transform (LT) in which a quadrature rule is applied on a suitable complex contour, 
namely a parabola. Methods of this kind have been successfully applied [13, 40] to 
the ML function restricted to some very special cases (0<a<l,/3 = lor real z). 

The extension to the more general case is however not trivial and demands not 
only a different and more in-depth theoretical analysis but also a thorough different 
strategy. Since the possible presence of a large number of singularities of the LT, 
it can indeed be impossible to find a contour encompassing all the singularities and 
behaving in a satisfactory way for computational purposes. Our approach is therefore 
to consider separate regions in which the LT is analytic and look, in each region, 
for the contour and discretization parameters allowing to achieve a given accuracy. 
The optimal parabolic contour (OPC) algorithm hence selects the region in which the 
numerical inversion of the LT is actually performed by choosing the one in which both 
the computational effort and the errors are minimized. 

The paper is organized as follows. Section 2 introduces the LT of the ML function, 
describes its analyticity properties and discusses the numerical inversion. In Section 
3 the OPC algorithm is presented and a detailed error analysis is derived in order 
to provide information for the selection of the optimal contour and of the suitable 
quadrature parameters. Section 4 is hence devoted to illustrate numerical experiments 
and some concluding remarks are discussed in Section 5. 

2. Evaluation of the ML function by LT inversion. During the last decades, 
an increasing amount of attention has been devoted to methods for computing special 
functions by inverting the LT; there are some key factors accounting for this interest: 

1. for several functions (including the ML) the LT has an analytical formulation 
which is much more simple than the function itself; 
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2. algorithms for the numerical inversion of the LT are usually quite simple to 
implement and run in a fast way; 

3. it is possible to derive accurate error estimations and perform the computation 
virtually within any prescribed accuracy. 

Although from a theoretical point of view the inversion of the LT is an ill-posed 
problem, satisfactory numerical results are expected for the ML function since it is 
possible to evaluate its LT in the whole complex plane and with high accuracy. 

An explicit representation of the LT of (1.1) and (1.2) is, however, not available. 
We must therefore introduce the following generalization of the ML function (1.1) 

(2.1) AeC, 
in order to express the corresponding LT as [20, 22, 29] 

£a,0{s]\) = < 1 ) 

(for easy of presentation we just focus on the two parameter function (1.1); the ex¬ 
tension to three parameter case (1.2) will be discussed in the Subsection 3.4). 

By means of the formula for the inversion of the LT it is possible to formulate 
the following integral representation of A) 

/*tT + ZOO 

( 2 . 2 ) ea,0{t;X) = — e’'^£a,p{s]\)ds, 

Z'K'l J(j—ioQ 

where (cr — ioo, a ioo) is the Bromwich line, with cr £ R chosen to ensure that all 
the singularities of £a,p{s\ A) lie to the left of the line 5J(s) = a. Since the presence of 
non integer powers, £’q_^(s; A) is a multi-valued function and a branch-cut extending 
from 0 to —oo along the real axis is introduced to make the integrand single-valued. 

Remark 2.1. For convenience we assume A ^ 0; it is readily verified that 
ea./3(^;0) = t^“^/r(/3). Moreover, also t = 0 is of no interest since ea,/3(0; A) = 0 for 
fi > 1, eap(0; A) = l/r(/3) and ea^pifi A) —?► -l-oo as t —>■ 0+ for fi < 1. 

As first suggested by Talbot [34], to exploit (2.2) for numerical computation it is 
necessary to deform the Bromwich line into an equivalent contour C that begins and 
ends in the left half of the complex plane in order to rapidly dampen the exponential 
factor e'** and avoid high oscillations which are source of numerical instability (for the 
equivalence of the contours it is meant that they encompass all the singularities of 
£a,p{s]\) to the left). Once a suitable contour is selected, a quadrature rule can be 
applied. 

The above two steps are intimately related. As deeply studied in [35, 40], the 
choice of the contour affects in a significant way the convergence properties of the 
quadrature rule which depend on the analyticity of the integrand in a region sur¬ 
rounding the path of integration. A satisfactory selection of the deformed contour is 
therefore not possible without a subtle analysis of the regions in which £a,p{s;X) is 
analytic. 

After denoting 9 = Arg(A), —tt < 9 < tt, the poles of £a^p{s-, A), i.e. the solutions 
of the equation s“ — A = 0, are 

(2.3) s* = Ai/“ = |A|i/“e*^, j £ Z. 
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The relevant poles are those in the main Riemann sheet, for which it is —tt < 
{9 + 2jTT)la < TT or, equivalently, such that j belongs to 


(2.4) 


J{a,e) = Ij €1. 


a 9 . a 9 

~ 2 ~ ^ - 2 ~ ^ 


their number depends on a and 0, ranging from zero when 0 < a < 1 and |0| > ai: 
to a possible very large number otherwise. 

The origin is a pole only when P > a; however, it must be always included among 
the singularities because a branch-point singularity occurs at the origin. 

From a formal point of view we denote with S* = {sg, s*,..., Sj} the set of 
all singularities (the poles and the branch-point) of £a,p{s;X), where sj = 0 and 
si = Jic., 9 y j = h with J = IJ(a,0)|. 

In the presence of a large number of singularities, or when some of them have 
large imaginary part, it can result nearly impossible to find suitable contours allowing 
a fast decay of the exponential factor and, at the same time, encompassing all the 
singularities. For this reason it can be useful, thanks to the Cauchy’s residue theorem, 
to remove some of the poles by residue subtraction 


(2.5) 


ea,p{t]\)= ^ Res(e®‘fa,/ 3 (s; A),S*) 



e’^*£a,i3{s;X)ds, 


where C S* is the set of all singularities of £a ,/3 laying on the rightmost part of 
the complex plane delimited by C and Res[f, s*) denotes the residue of the function 
/ at s* (observe that, due to the selected branch-cut, C cannot traverse the negative 
real semi axis and must encompass at least Sg = 0 to its left). 

It is a favorable achievement that the residues in (2.5) can be explicitly represented 
in terms of elementary functions as 

i?es(e**£:„,;3(s;A),s*) = 

Assumed the contour C be represented by means of a complex-valued function 
z(u), —oo < u < oo, then (2.5) can be rewritten as 


( 2 . 6 ) 


^a,/3 (j'l '^) — 

a 


(^*) 


s* ^Sp> 


1 


where 

(2.7) 


g{u) = e^(“)‘fa./ 3 (z(M); A)z'(m) 


(z(m))“ - A 


Numerical quadratures for integrals on unbounded intervals 



are presented in several papers and reference books (e.g., see [6]). An extensive 
analysis of the trapezoidal rule has been recently provided in the remarkable paper 
by Trefethen and Weideman [35] which not only focuses on the fast convergence of 
the trapezoidal rule but also discusses its main practical applications. Despite its 
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simplicity, the trapezoidal rule appears indeed as a powerful tool to perform fast and 
highly accurate integration in a variety of applications. 

On a given equispaced grid kh, k G 'Z, with step-size h > 0, the infinite and finite 
trapezoidal approximations of / are 

OO N 

Ih = h g{kh), Ih,N = ^ X] 9 {kh) 

k— — oo k— — N 

and the corresponding error I — results from the sum of the discretization error 
DE = |j — /;i| and the truncation error TE = \lh — Ih,N\- 

Under the assumption that g(u) decays rapidly as u —> ±oo, an estimation of TE 
is given by the last term retained in the summation, i.e. TE = 0(\g{hN)\^, N —>■ oo. 

As discussed in [35, 40], the estimation of DE is performed on the basis of the 
analyticity properties of g(u). For reasons which will be clear later, we need to slightly 
modify the statement of the original result, with no substantial changes in the proof 
which remains the same as outlined in [35]. 

Theorem 2.2. Let g{w) be analytic in a strip —d* < '^(w) < c*, for some c* > 0, 
d* > 0, with g{w) —>■ 0 uniformly as jwj oo in that strip. For any 0 < c < c* and 
0 < d < d* it is 


where 


and 


DE=\I-Ih \ < DE+{c) + DE_{d), 




/ oo pCO 

\q(u + ir)\du, M-(d) = max / \q(u — ir)\du. 

-oo' ' 0<r<dJ_J^^ 


In most cases (for instance with the exponential function [40]), the contribution of 
M+(c) and M-{d) is negligible and the estimations DEj^ic) « DE+{d) ~ 

g-27rd/?i sufficiently accurate for a satisfactory error analysis. When applied to the 
ML function it is possible, depending on the parameters a and /3, that M+(c) —>■ +oo 
when c —?► c* and M-{d) -boo when d ^ d*. The consequence of unbounded 
limits for M+{c) and M-{d) is that their contribution can be nonnegligible. This 
is especially true within very narrow strips of analyticity (as when there are several 
singularities), for which c or d are necessarily close to their upper bounds c* and d*. 

Providing a reliable estimation for M+ (c) and M_ (d) (and for the rate by which 
they tend to +oo) and including them in the error analysis is therefore of utmost 
importance in order to select optimal parameters and fulfill an assigned accuracy. 

3. Parabolic contours and the OPC algorithm. Removing some of the 
singularities by the residue subtraction in (2.5) offers a considerable freedom in the 
choice of the integration path. The task of selecting a suitable contour in a specific 
region of the complex plane is greatly simplified by first fixing the geometric shape 
and hence adopting a parametrized description of the contour with very few (usually 
just one) parameters; the problem is thus reduced to the evaluation of the optimal 
parameters. 
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Several families of contours have been proposed so far. After the original work of 
Talbot [34] on contours of cotangent shape (see also [7, 28, 38]), a special attention 
has been paid to parabolic [2, 14, 39, 40] and hyperbolic contours [21, 14, 33, 40]. 

The convergence rates of the A^-points trapezoidal rule on cotangent, hyperbolic 
and parabolic contours have been studied in [36]; the respective rates of 0(3.89“^), 
0(3.20“^) and 0(2.85“^) indicate a fast convergence with all these contours. 

Although the convergence with cotangent and hyperbolic contours is slightly 
faster, the simpler representation of parabolic contours makes them much more easy 
to handle; therefore, parabolas appear to be preferable especially when the presence 
of a certain number of singularities demands the fulfillment of tightened constraints. 

As in [35, 40], for a real parameter /i > 0 we consider the parabolic contour 

(3.1) C : z(u) = fjiiyiu + l)^, —oo < u < oo. 

To select the singularities that must be removed in (2.5), we partition the com¬ 
plex plane in neighboring regions having the singularities of £a ,/3 on their respective 
boundaries; in each region the parabolic contour and the discretization parameters 
are determined, according to a suitably modified version of the procedure described 
in [13, 40], with the aim of fulfilling a prescribed accuracy £ > 0. Among the possible 
contours (one for each region), the OPC algorithm makes an optimal selection with 
respect to the computational effort: the region and the contour involving the smaller 
number N of quadrature nodes is selected. Nevertheless, some issues related to reduce 
the propagation of round-off errors are also addressed. 

As already observed in [13], the computation necessary to select the contour is 
much less than the computation involved by the inversion of the LT. Thus the overall 
process of establishing the contour in an optimal way adds only a negligible amount 
of computation, with the obvious advantage of performing the actual, and more ex¬ 
pensive, inversion with the smallest possible number of floating point operations. 

The starting step in the OPC algorithm is to sort the singularities of £a,p in order 
to identify a sequence of regions delimited by two consecutive singularities. To this 
purpose we introduce the function : C —t K+ defined according to 

The function (p allows to split the complex plane in regions bounded by parabolas 
of type (3.1) as stated in the following Proposition. 

Proposition 3.1. Let z(u) = fx(iu -l- 1)^, with p, > 0. A point s € C lies on the 
parabola z{u) whenever ip{s) = /r. Moreover, whenever (p{s) < p the point s lies at 
the left of the parabola z(u) and whenever ip{s) > p the point s lies at the right of the 
parabola z(u). 

Proof. After expanding z{u) = p{l — vf) + 2iup, it is immediate to see that a 
point s G C lies on the parabola described by z{u) whenever 

f p{l — vf) = 5i(s) 

2pu = 3(s) ’ 

hence, the first part of the proof follows by considering the unique positive solution p 
of this system. The remaining statements easily follow by observing that z{u) defines, 
as p varies, a bundle of parabolas moving from left towards right as p increases. □ 
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Hereafter, we assume that the singularities s* G S* are arranged in ascending 
order with respect to the function ip, i.e. 

0 = ip{sq) < v7(4) < ■ • ■ < ‘fiis}). 

We can therefore consider the J + 1 parabolas defined by (3.1), each with p, = 
ip{s*), and determining J + 1 disjointed regions in the complex plane with the 
singularities s* and on the left and on the right boundary (except for the last 
region Rj which is unbounded to the right). In Figure 1 we show, for instance, the 
complex plane partitioned into 6 regions of this kind (note that the first parabola, the 
one with fj, = ipis^) = 0, collapses onto the negative real axis). 



Fig. 1. Partitioning of €, into some regions Rj by means of parabolas (3.1) through s'). 


A parabolic contour defined according to (3.1) is in a region Rj, j G {0,1,..., J}, 
whenever p, satisfies 


(3.2) if{s*) < n< ip{s*+i). 

(for convenience, a fictitious singularity with (/^(sj^^) = +oo is introduced). At 
the same time, if Zj{u) = Pjiiu + 1)^ is a parabolic contour with pj satisfying (3.2), 
theni?j is the region of analyticity of Theorem 2.2 for (it) = e^^^^'>*£a,p{zj(u)-,\)z'j{u) 
and can be expressed as 


(3.3) = {z G C I z = fj,j(iw + 1)^, w G C and — d* < S(ii;) < c*} , 


with 


(3.4) c* = 


I 



J = 0 

i> 1 



j<J 
3 = J- 


Basically, by means of (3.4) it is possible to represent each region of analyticity 
Rj as the strip —d* < ^{w) < c* in the ic-plane. 

The estimation of the discretization error of Theorem 2.2 involves finding upper 
bounds for M_|_(cj) and M-{dj), with Cj < c* and dj < d*, in each region Rj. 
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Matching the corresponding DE^{cj) and DE_{dj) with TE = = 

0[ ), according to a procedure similar to that devised in [40], allows to 

obtain the optimal parameters (contour geometry fXj , step-size hj and number Nj of 
quadrature nodes) in order to achieve a prescribed tolerance e > 0. 

The region involving the minimum computational effort (i.e., the one with the 
minimum number of quadrature nodes) is hence selected to perform the numerical 
inversion of the LT; the residues corresponding to the singularities left out by the 
selected contour are accordingly added in the final result as stated by (2.5). 

The main steps of the OPC algorithm can be therefore listed as follows: 

1 . estimation of M+{cj) and M-{dj) in each region Rj; 

2. matching, in each region Rj, of DE+{cj), DE-{cj) and TE with the pre¬ 
scribed accuracy e > 0 and evaluation of the parameters , hj and Nj ; 

3. selection of the region Rj in which to perform the integration on the basis of 
the lowest computation and reduction of round-off errors. 


3.1. Estimation of Mj^{cj) and M_{dj) in each region Rj. To provide an 
estimation of Mj.(cj) we distinguish the case in which the region Rj is bounded to the 
left by the singularity at the origin (i.e., j = 0) and the case in which the singularity 
on the left boundary of Rj is one of the poles of Ea,f 3 except the origin (i.e., j >0). 

To discuss the first case we introduce the following preliminary result. 

Lemma 3.2. Let A,a > 0 and p G R. Then as A —>■ 0 



'{u^ + Af du = 


oil) 

O(log a A) 
0{AP+^) 


ifP > 
ifP=-l 
ifP < - 5 - 


Proof. By splitting the integral into the two subintervals (—oo, 0] and [0, oo) and 
making the change of variable s = u^/A, it is possible to preliminarily observe that 



{u^ + Af du = AP+i 



s-5e-‘"^"(s-f l)^ds. 


The right-hand side of the above equation is the integral representation of the con¬ 
fluent hypergeometric function of the second kind [37], namely the 'I'(a, b, z) function 
with parameters a = 1/2, 5 = p-|-3/2 and z = crA, and hence 



-I- A)^ du 


AP+5r(i)vi/(i,p + |,aA). 


As z —?► 0 the T-function admits the following asymptotic expansions [42] 

^ OO OO 

if6^Z 


4'(a, 6, z) = <' + z 


j=0 

oo 

3=0 j=0 


(1) 7 
v) >Z^ 


3=0 

-b 


(2) 3 
v) 


+ log(z)zi-^ if 6 e Z- U {0} 


h-1 


^ uf z^-+ ^ uf z-^-+ log(z) ^ 


3=0 


(3) 3 
Wj z-' 


if 6 G Z+ 


3=0 


1=1 


1=0 


with coefficients {u^P}j, {Vj^^}j and {w^^^}j, i = 1,2,3, independent of z. The proof 
now follows by considering the leading terms in each summation. □ 
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The first region Rq is bounded to the left by the singularity Sq at the origin; since 
(3.4), the corresponding upper bound for the strip of analyticity in the w-plane is 
Cq = 1. We can provide the following estimation for M+(co) for cq < 1. 

Proposition 3.3. Let fj,o be sueh that 0 < /iq < f'or any cq < 1 there 

exists a constant M_|_ (independent of cq) such that 

M+(co) < M+ • M(co), 

where as 1 it is 

r oil) ifP<a + l 

M{c) = l C>(log/xot(l - c)2) = e>(log(l - c)) i//3 = a + l 

[ 0((l-c)2(“-/3+i)) f//3>a + l. 


Proof. By replacing (3.1) and z'iu) = 2/ro(j — u) in (2.7), we preliminarily obtain 

^-P+l ..ua((l-r)+zuPt ((1 - ^) + + ^ 


g{u + ir) = 2ifiQ' 


z{u + ir)°‘ — A 


Since A 0 (see Remark 2.1), it is natural to assume the existence of a positive 
M such that \z{u + ir)°‘ — X\ > M for any r G [0, Cq). Hence 


\giu + ir)| < {u^ + (1 - 

and, after putting for shortness M+ = we have for any cq < 1 

/ OO 

e“^o“ t (y 2 -I- (1 _ r)^)“ ^ du. 

-OO 


The proof now follows after applying Lemma 3.2. □ 

With A very close to 0, in the above proof it is possible that M <C 1, thus affecting 
the asymptotic estimation for M(c). In this case, we are in the presence of a very 
narrow region Rq which, as we will discuss later in the final part of subsection 3.2.1, 
must be discarded since it does not allow one to achieve the assigned tolerance. For 
this reason, we do not consider the effects on M{c) of a possibly very small M. 

We now consider the regions Rj, j = 1,... ,J, which are bounded to the left by 
one of the poles s* of £a,p{s', A) except the origin (i.e., <p{s*) > 0). 

Lemma 3.4. Let a, 6 G K, with b > a > 0. For 0 < a < 1 it is 6“ — a“ > 
ab°‘~^{b — a) and for a > 1 it is b'^ — 0 °^ > aa°‘~^ib — a). 

Proof. For 0 < a < 1 it is immediate to verify that 


6“ - a“ = a 


„Q!— 1 


ds > a 6“ ^ ds = ab^^ ^{b — a) 


and in a similar way the proof follows for a > 1. □ 

Proposition 3.5. Let j g {1,..., J}, gj > 0 such that pis*) < gj < <p(s*_|_]^) 
and c* the upper bound of the strip of analyticity (3.3) corresponding to Rj. For any 
Cj < c* there exists M+ > 0 (independent ofcj) such that 


M+ic,) < M+ ■ [c* - c,)-\ 
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Proof. Let Cj < c* and consider r G [0, Cj]. By Proposition 3.1, s* lies on the 
parabola z*(u) = (p(s*){iu + 1)^ and hence A = (sp“ lies on the curve = 

{iu + 1 )^“. The distance between (zj(u + ir))°‘ and A is therefore greater 
than the distance between {zj(u + ir))“ and ( 2 :*(m))“ evaluated at u = 0 , i.e. 


\{zj{u + ir))^ - X\ > 


Since <p{s*) < — r)^, Lemma 3.4 yields 

\{zj{u + ir))“ - A| > aPaj (l - , 

where Paj ^ when 0 < a < 1 and Paj = /r““^(l — r)^“ ^ when a > 1 . 

It is elementary to see that for any r < Cj < c* it is 

(^(1 - = ^^o ((l - - (1 - cfjf) 

= - r){2 - r - c*) > 2fij{c* - Cj){l - c*) 

Since form (3.4) {p{s*)/^j = (1 — c*)^, when 0 < a < 1 we can easily verify that 

\{zj{u + ir))°‘-\\ > 2a((/j(s*))“"Vi(Cj-Cj)(l-cp = 2a{ip{sy)°‘ {c* - Cj){l- c*y~^ , 

while, for a > 1 , it is instead 

\zj{u + irY - A| > 2a^“(c* - Cj){l - 

We denote with Qa,j the constant independent of Cj 

2a((^(s*))“(l-c*)-i if0<a<l 
\ 2a/Lt“(l - if a > 1 

and write the inequality 

\{zj{u + ir))°‘ - A| > Qa,j{c* - Cj) 

which allows to obtain 

noo 

M+{cj) = max / \g[u + ir)\du < M+{c* — Cj)~^, 


with M+ = and 

M+= —— max^M(r), M{r) = / e“'^“ * (m^ + (1 — du 


being M(r) < +oo since r < c* < 1 . □ 

Also for M-{dj) it is necessary to distinguish two cases: when the computation 
is performed in a region i?^, j = 0,.. ., J — 1, bounded to the right and when the 
integration is instead performed in last region Rj. 

Proposition 3.6. Let j € {0,..., J — 1 }, fij > 0 such that f{s*) < g,j < 
and d* the lower bound of the strip of analyticity (3.3) corresponding to Rj. For any 
dj < dj there exists M_ > 0 independent of dj) such that 

M.{d,) < M_ . (d* - d,)-\ 
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Proof. The prof is symmetric to the proof of Proposition 3.5 and we omit the de¬ 
tails. We just point out that the term M_ is now given by M_ = 
with M_ obtained in a similar way as M+. □ 

The discussion for M_(dj) in the last (right-unbounded) region Rj is the same 
as proposed in [40]. We just recall that an upper bound for M-{dj) is achieved at 
dj = Tr/ifijthj) — 1 which allows to write 

(3.5) i:)£;_(dj) = hj^O. 

3.2. Evaluation of the quadrature parameters. Thanks to the analysis car¬ 
ried out in Subsection 3.1, and after highlighting the exponential growing factor in 
M -, some upper bounds for the discretization errors are now available in the form 

DE+{cj) < DE_{dj) < 

for Cj < c*, dj < d* and some nonnegative values pj and only in the last region 
Rj a different result applies for DE-{dj), according to (3.5). 

Unless Pj = 0 and qj = 0, the presence of the algebraic terms (c* — Cj)~P^ and 
{dj — dj)~'^^ cannnot be disregarded; they can be indeed very large and an unfit 
selection of Cj and dj can lead to an incorrect error analysis as already observed in 
the first part of [13]. 

The task of including, in the error analysis, the contribution of a possibly large 
algebraic term was accomplished in [13] by introducing an auxiliary variable and 
expressing the parameters in the formula for the numerical inversion of the LT in 
terms of this variable; the optimal value of the auxiliary variable was hence selected 
by minimizing the number of quadrature nodes in order to keep the computational 
effort at a minimum. Numerical experiments showed that despite the nonnegligible 
computation required by finding the minimum of a nonlinear function, the overall 
computation was the same performed in an efficient way. 

The work in [13] was anyway devoted to the evaluation of the ML function (2.1) 
on the real negative semiaxis and, mainly, for 0 < a < 1, thus involving just one 
singularity, namely at origin. In the more general context this approach is no more 
feasible: since most of the regions Rj have two distinct singularities, the introduction 
of two auxiliary variables leads to the need of finding the minimum of a nonlinear 
function with two variables, a problem whose solution can be quite expensive. 

We propose here a completely different approach to take into account the algebraic 
terms in DE^{cj) and DE-{dj). To this purpose we distinguish again two main cases: 
the case of a bounded to the right region Rj (i.e., j = 0,..., J — 1) and the case of 
the right-unbounded region (i.e., the last region Rj). 

3.2.1. Quadrature parameters in a region bounded to the right. The 

most straightforward way to prevent the possible growth of the terms (cJ — Cj)~P^ 
and {d* — dj)~‘^^ in the discretization errors DE+{cj) and DE-{dj) is by forcing pj 
to belong to a subinterval [(p*j,(p*j_^-^\ instead of [i,5(s*), :/9(s*_|_j)], with 

(3.6) ip{s*j) <(p*j < (p*j+i < p{s*+i) 

(the equality in (3.6) is introduced just to cover the case in which pj =0). Under the 
conformal map Zj{w) = pj{iw + 1)^, the strip —d* < 3(w) < c* corresponding to the 
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region determined by the parabolas through ip* and ifj^i is given by 


(3.7) 







A suitable strategy is to fix an arbitrary value / > 1 and select ip* and ‘P*j by 
forcing 

(3.8) (c* - -c*)-^^ = (d* - d*)-^^ = f 


Whenever (3.8) can be imposed without a large value of /, for instance / ~ 1, 
the terms (cj — Cj)~Pi> and {d* — d*)~‘^^ can be neglected in DE+{cj) and DE-{dj). 
Within very narrow regions Rj it cannot be possible to satisfy (3.8) for sufficiently 
small values of /; the accuracy e must be therefore scaled as e = e// before removing 
{c* — c*)“P^' and {d* — d*)~^i> from the discretization errors. Obviously, this scaling 
also has an effect on the truncation error TE and leads to slightly increase the number 
Nj of quadrature nodes. 

The asymptotic balancing of the different components of the error now reads 


(3.9) 


hj 


^TTd* 


I^Nf) 


loge 


from which it is immediate to obtain 


(3.10) 


2 + w j 


w = — 


^ 1 + 1 ^ 

loge ’ 


and 


27r 

loge (l + w)y^ + 





loge 
td-3 ' 


An essential task is to select / small enough to make {c* — c*)~Pii and {d* — d*)~'^^ 
negligible in DE+{cj) and DE-{dj) and, at the same time, satisfy (3.6). To this 
purpose we explicitly represent c* — c* and d* — d* in terms of (p* and as 

* (2 + ^) (y/^-(2 + n;)(^^(s*+i)-^^) 

(l + ^)y^+y^ ’ J J- + ■ 


The obvious assumption / > 1 is sufficient to ensure that p^{sj) < (p*j and 'P*j+i < 
‘/5(s*_|_i); anyway, a minimum threshold value /mm > 1 must be determined with the 
aim of fulfilling (p* < ip*j^-^ for / > /mm. 

When pj = 0, a simple computation allows to verify that ip{s*) = ip* < for 


fmin — 
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in the more general case we provide the following result (note that in regions bounded 
to the right it is always qj 7 ^ 0 ). 

Proposition 3.7. Let Pj,qj > 0 and rj = ma.x{pj, qj} . If 


(3.11) 


f>fn 


fmin — 




then (p* < plj+i ■ 

Proof. It is elementary to verify that (c* — = {d* — d*)~'^^ = f when (p* 

and are obtained after solving the linear system 


2 + w- {l+w)f-^/P^ -f-^IPi 
(l + w)/-!/-?^ 2PwPf-^l^- 
whose solutions can be explicitly formulated as 




= {2+w) 


= 


{2 + W + ) ^(p{s*) + f~^IP^ 




/ 




and 


P*j+i 


2 + ’w—{l + w)f ^/Pi + f l/lj 

-(1 + ^p{s*) + {2 + w-{l + 

2 + W — (1 + w)f~^^P3 + 


By the hypothesis (3.11) it is 


/ 


-I/Pi 


< 


'i+1 


< 1 , / 


-i/ij 


< 


^l^) + ^fp{^) 


and hence 

(3.12) /-i/«^Y(^(s*) + /-i/«y^(^(s*+i) < - ^<^(s*+i)- 

A simple computation allows to prove that 


< 1 


7’i — P*j+i + (2 + w) 


2 + W — (1 + w)f~^/Pi + 


from which the proof follows after using (3.12). □ 

In very narrow regions the value of p{s^) can be very close to p{s j+i) and the 
threshold fmin can be too large to assure the achievement of a small tolerance e > 0 ; 
in this case no contours can be selected and the region must be discarded. 

Remark 3.8. A more eonservative error analysis would take into account also 
the exponential growing term e^A in M_|_ (see the proof of Propositions 3.3 and 3.5). 
In this case, and by using for simplicity the upper hound e'^^+A, the integration pa¬ 
rameters are obtained after replacing loge with loge — ip)+it in the formulas for w 
and hj. This change however does not seem to offer substantial improvements sinee it 
aetually exerts the effects in regions with large pis*) and <p(s*_|_]^) whieh are normally 
discarded for aeeuraey reasons as we will diseuss later on. 
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3.2.2. Quadrature parameters in an unbounded region to the right. In 

the last and right-unbounded region Rj the balancing of the exponential factors of 
the errors leads to 


(3.13) 


27rG) 


27r 


+ — = - h^jNj) = loge, 


hj hj 

where Cj < Cj is selected according to (3.7) for (p*j > ip(sj), from which we obtain 


1 _l_ 

(3.14) hj = f,j = 

tv,7 


ttNj 


2t(l + c})(l + 2^j) 


, iVj = -i±^loge. 


Unfortunately, because of the implicit dependence on the unknown /ij, we cannot 
use (3.7) to determine Cj. It is therefore necessary to formulate hj, fxj and Nj directly 
in terms of (p*j instead of Since it is 

(3.15) ,, = 

by matching the two representations of fij in (3.14) and (3.15) we obtain a second 
order algebraic equation with respect to Cj whose unique solution satisfying Cj < 1 is 

3 + ^ - VTTTTA ^ ttNj 
'3 =- ^ = 

A straightforward manipulation leads to 

^ 1 / 3A 2-2Vl + 12A\ 

V 4-A 4-A y ’ (7-Vl + 12ij^ 

and, after imposing —2TTc*jlhj = loge, we obtain 


Nj 


l + 2 cy 
27rcy 


loge 



3logs , I, plogg^ 

2^*jt Y t(p*j )' 


A direct evaluation of a suitable value for (p*j is now not possible since its implicit 
dependence on /ij. We hence use an iterative process by which, starting from an 
initial guess very close to Lp{s*j), the value of (p*j is increased until the corresponding 
value of /, evaluated as 


(3.16) 


f — ~ Cj) 


\ ) 


does not fall into an interval [fmin, fmax] which is a priori selected, for instance 
[1,10]. To this aim a target value ftar G [fmin, fmax] can be established and the 
new attempted value for (pj is obtained after replacing / with ftar in (3.16). In our 
experiments we have observed the convergence of this procedure in very few (usually 
1 or 2) iterations. 
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3.3. Selection of the region in which to invert the LT. After evaluating 
parameters fij, hj and Nj in each subregion Rj, we select the region involving the 
minimum number Nj of quadrature nodes to actually perform the numerical inversion 
of the LT with the minimum computational effort. 

Because of the presence of the factor , with large values of t and/or ^ it 
is possible the presence in the summation Ih,N of terms with large magnitude and 
terms with small magnitude; the effects of this simultaneous presence are in numerical 
cancellation which can become catastrophic. 

As already observed in [39], the rounding error is roughly RE Ri e^‘e, with e the 
machine precision. To keep rounding errors below the desired accuracy e > e it is 
therefore necessary that ^j < (loge — loge)/t and, hence, from (3.10) it is sufficient 
to verify 



In regions with 'f{Sj) > (log£ — loge)/t this condition cannot be fulfilled; in order 
prevent round-off errors from destroying all of the significant digits, such regions must 
be discarded and the computation moved to one of the remaining regions. In the other 
cases the above equation provides a bound for (p*j+i- 

Another possible source for numerical cancellation is the closeness of the contour 
to one of the singularities on the boundary of the region Rj. We observe however that, 
despite the previous case in which the accuracy is affected by a factor proportional to 
e^*e, in this case the accuracy is affected only in an algebraic way and, as observed 
by means of numerical experiments, it is sufficient to select (p* and ‘P*jj.i as previously 
described in order to avoid the cancellation. 

In the last region Rj it is possible, even when </?(sj) > (log£ — loge)/t, that the 
value pLj resulting from the balancing of the various components of the error is too 
large and the round-off error RE k. e^-'^e exceeds the required tolerance e. Since in 
this case RE dominates the discretization error DE- [39], it is necessary to replace 
in (3.13) the exponential factor of DE_ with that of RE] by solving explicitly with 
respect to fij, hj and Nj we derive in this case 



The introduction of the round-off error in the error analysis prevents from placing 
the contour in a place in which it is not possible to guarantee that round-off errors 
do not exceed the assigned tolerance e. Since now an explicit value of /rj is available, 
the computation of (p*j follows directly from (3.16) as 



(obviously, the region Rj must be discarded when (pj exceeds the threshold (log e — 
loge)/t since even the accuracy fe cannot be achieved). 

Since rounding errors depend, in an exponential way, on the value of t, it can be 
useful to scale the ML function in order to force t to assume small values. By simple 
algebraic manipulations, it is indeed immediate to see that for any r > 0 



and it is therefore possible to reduce the propagation of rounding errors by suitably 
using the above scaling, for instance with t k. t. 
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3.4. Extension to three parameters. The main information used by the OPC 
method is the location and the strength of the singularities of the LT; its extension 
to the 3 parameter ML function (1.2) is therefore straightforward. The LT of the 
corresponding generalization A) = is indeed 

„a7-/3 

^ 2 ./ 3 (s; _ ;^)7 ’ ^(■5) > 0 and |As““| < 1, 

which has the same singularities of the 2 parameter counterpart. It is elementary to 
reformulate Proposition 3.3 by replacing a with ay and evaluate the new bounds of 
Propositions 3.5 and 3.6 respectively as M+(cj) < M+ ■ (c* — Cj)~'^ and M-{dj) < 
M_ • (d* - d,)-E 

With 7 1 we must restrict the computation to 0 < a < 1 and | Arg(A)| > air 

since otherwise non trivial difficulties (whose discussion is beyond the scope of the 
present paper) arise due to more involved branch-cuts; the case 0 < a < 1 and A real 
and negative is however the most interesting for applications [3]. 

4. Numerical experiments. To test the proposed method and verify its com¬ 
putational efficiency we present in this Section some numerical experiments. 

All the experiments are performed in Matlab, version 7.9.0.529, on an Intel Dual 
Core E5400 processor running at 2.70 GHz under the Windows XP operating system; 
the Matlab code implementing the OPC method and described in the previous sections 
is made available at [10]. As reference we use the values evaluated after summing the 
series (1.1) or (1.2) in variable precision arithmetic with 100 digits by means of Maple. 

In all the experiments we set the target tolerance e = 10“^®; the goal is to test 
whether it is possible to provide an approximation E'^ p{z) of the ML function ^{z) 
with an accuracy very close to the machine precision. The tolerance e represents the 
absolute error in the computation of the integral in ( 2 . 6 ) and this is the error we 
expect when the value of the function is not large in modulus (in this case no residue 
calculation is usually involved); otherwise, the summation of residues can dominate 
the integral in ( 2 . 6 ) by several orders of magnitude and the leading error is the round¬ 
off error in the computation of residues: in the double precision used by Matlab it 
involves a relative error smaller than e = 10“^®. The resulting error is therefore a 
combination of absolute (with small values of E^ p{z)) and relative (for large values 
of E2 errors and it can be represented as 


(4.1) 


- El,{z)\ 

^ + KA^))\ - 


In Figure 2 we report the error (4.1) for the 2 parameter function Ea^pi^z)^ for 
a = 0.7 and /? = 1.0, evaluated in several points z on the real negative axis. As we 
can clearly see, the OPC method achieves an accuracy very close to or smaller than 
the requested tolerance of 10“^^ (the few gaps in the error plot are due to the fact 
that in some cases the approximated and reference values are exactly the same). 

To show the efficiency of the proposed method we present in Figure 3 the com¬ 
putational time and we compare it with that of the Matlab mlf code [30]. This is so 
far the unique available Matlab code for the ML function and, since it is widely used, 
it can be considered as a sort of benchmark for testing new methods. 

We observe that whilst the CPU time consumed by OPC remains nearly constant, 
the mlf code demands for a CPU time close or slightly less than OPC for very small 
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Fig. 2. Error for E^^piz) with a = 0.7, /3 = 1.0 and arg(z) = tt. 


and large values ot \z\ whilst for moderate values of | 2 ;| the CPU time of mlf is some 
order of magnitude higher than OPC. 



— OPC 

— -mlf 


20 


Fig. 3. Computation time for Ea^^{z) with a = 0.7, /3 = 1.0 and £LTg(z) = tt. 


This nonuniform behavior can be explained by observing that mlf uses different 
techniques according to the value of \z\: for very small \z\ the series ( 1 . 1 ) is evalu¬ 
ated until numerical convergence and this computation is quite fast; an asymptotic 
expansion is instead used when |z| is large and the computation becomes faster and 
faster as |z| grows; for intermediate values of \z\ a Romberg integration is applied to 
an integral representation of the ML function, with a computational cost proportional 
to 2P whenever an accuracy e = 10“^ is requested. On the other hand, most of the 
computation of OPC is spent by the trapezoidal rule whose cost depends essentially 
on the number of nodes which is kept at the minimum by the algorithm (and it is 
roughly proportional to p for any argument z); the amount of computation required 
by the other tasks of OPC, such as location of the singularities, choice of the suitable 
region and evaluation of the quadrature parameters, is usually negligible. 

The plot in Figure 4 shows that the OPC algorithm behaves in a robust way 
and provides results within the requested tolerance also for complex values on the 
imaginary axis (we used here a = 0.5, j3 = 1.0 for which it is known that mlf does 
not provide accurate results). 

We conclude our experiments by presenting the errors for the three parameter 
function E'^ p{z) for a = 0.6, fi = 0.9, 7 = 1.2 and arg(z) = As we can see from 
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Fig. 4. Error for E(y^p[z) with a = 0.5, /3 = 1.0 and arg( 2 :) = 


Figure 5, OPC behaves in a satisfactory way and produces errors very close to the 
target tolerance also for E'^ ^{z). 



Fig. 5. Error for p{z) with a = 0.6, /3 = 0.9, 7 = 1.2 and arg( 2 :) = 

We do not report the CPU time for E'^ p {z) since it would not provide any further 
information; as discussed in the Subsection 3.4, the evaluation of the three parameter 
function just involves different coefficients in the error estimations and most of the 
computation (and hence the CPU time) is the same as in the two parameter case. 

5. Concluding remarks. In this work we have presented the OPC method for 
the evaluation of the two parameter ML function, a function which plays a fundamen¬ 
tal role in fractional calculus. The OPC method allows to evaluate the ML function 
with high accuracy and numerical experiments have shown its computational effi¬ 
ciency. The generalization to the three parameter ML function has been discussed 
and tested too. The corresponding Matlab code is made freely available [10]. 

Acknowledgments. The author is extremely grateful to the anonymous referees 
for their insightful and constructive remarks which allowed to improve the paper in a 
remarkable way. 
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